/* Sound.cpp
 *
 * Copyright (C) 1992-2024 Paul Boersma
 *
 * This code is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This code is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
 * See the GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this work. If not, see <http://www.gnu.org/licenses/>.
 */

/*
 * a selection of changes:
 * pb 2006/12/31 stereo
 * pb 2010/03/26 Sounds_convolve, Sounds_crossCorrelate, Sound_autocorrelate
 */

#include "Sound.h"
#include "Sound_extensions.h"
#include "NUM2.h"

#include "enums_getText.h"
#include "Sound_enums.h"
#include "enums_getValue.h"
#include "Sound_enums.h"

Thing_implement (Sound, Vector, 2);

void structSound :: v1_info () {
	structDaata :: v1_info ();
	MelderInfo_writeLine (U"Number of channels: ", our ny, our ny == 1 ? U" (mono)" : our ny == 2 ? U" (stereo)" : U"");
	MelderInfo_writeLine (U"Time domain:");
	MelderInfo_writeLine (U"   Start time: ", our xmin, U" seconds");
	MelderInfo_writeLine (U"   End time: ", our xmax, U" seconds");
	MelderInfo_writeLine (U"   Total duration: ", our xmax - our xmin, U" seconds");
	MelderInfo_writeLine (U"Time sampling:");
	MelderInfo_writeLine (U"   Number of samples: ", our nx);
	MelderInfo_writeLine (U"   Sampling period: ", our dx, U" seconds");
	MelderInfo_writeLine (U"   Sampling frequency: ", Melder_single (1.0 / our dx), U" Hz");
	MelderInfo_writeLine (U"   First sample centred at: ", our x1, U" seconds");
	integer numberOfCells = our nx * our ny;
	bool thereAreEnoughObservationsToComputeFirstOrderOverallStatistics = ( numberOfCells >= 1 );
	if (thereAreEnoughObservationsToComputeFirstOrderOverallStatistics) {
		double minimum_Pa = our z [1] [1], maximum_Pa = minimum_Pa;
		longdouble sum_Pa = 0.0, sumOfSquares_Pa2 = 0.0;
		for (integer channel = 1; channel <= our ny; channel ++) {
			constVEC const& waveform_Pa = our z [channel];
			for (integer i = 1; i <= our nx; i ++) {
				const double value_Pa = waveform_Pa [i];
				sum_Pa += value_Pa;
				sumOfSquares_Pa2 += value_Pa * value_Pa;
				if (value_Pa < minimum_Pa)
					minimum_Pa = value_Pa;
				if (value_Pa > maximum_Pa)
					maximum_Pa = value_Pa;
			}
		}
		MelderInfo_writeLine (U"Amplitude:");
		MelderInfo_writeLine (U"   Minimum: ", Melder_single (minimum_Pa), U" Pascal");
		MelderInfo_writeLine (U"   Maximum: ", Melder_single (maximum_Pa), U" Pascal");
		const double mean_Pa = (double) sum_Pa / numberOfCells;
		MelderInfo_writeLine (U"   Mean: ", Melder_single (mean_Pa), U" Pascal");
		const double meanSquare_Pa2 = (double) sumOfSquares_Pa2 / numberOfCells;
		const double rootMeanSquare_Pa = sqrt (meanSquare_Pa2);
		MelderInfo_writeLine (U"   Root-mean-square: ", Melder_single (rootMeanSquare_Pa), U" Pascal");
		const double energy_Pa2s = (double) sumOfSquares_Pa2 * our dx / our ny;   // Pa2 s = kg2 m-2 s-3
		MelderInfo_write (U"Total energy: ", Melder_single (energy_Pa2s), U" Pascal\u00B2 sec");
		const double rho_c = 400.0;   // rho = 1.14 kg m-3; c = 353 m s-1; [rho c] = kg m-2 s-1
		const double energy_J_m2 = energy_Pa2s / rho_c;   // kg s-2 = Joule m-2
		MelderInfo_writeLine (U" (energy in air: ", Melder_single (energy_J_m2), U" Joule/m\u00B2)");
		const double physicalDuration_s = our dx * our nx;
		const double power_W_m2 = energy_J_m2 / physicalDuration_s;   // kg s-3 = Watt/m2
		MelderInfo_write (U"Mean power (intensity) in air: ", Melder_single (power_W_m2), U" Watt/m\u00B2");
		if (power_W_m2 != 0.0) {
			const double referencePower_W_m2 = 1.0e-12;   // this equals the square of 2.0e-5 Pa, divided by rho c
			const double power_dB = 10.0 * log10 (power_W_m2 / referencePower_W_m2);
			MelderInfo_writeLine (U" = ", Melder_half (power_dB), U" dB");
		} else {
			MelderInfo_writeLine (U"");
		}
	}
	const bool thereAreEnoughObservationsToComputeSecondOrderChannelStatistics = ( our nx >= 2 );
	if (thereAreEnoughObservationsToComputeSecondOrderChannelStatistics) {
		for (integer channel = 1; channel <= our ny; channel ++) {
			const double stdev = NUMstdev (our z [channel]);
			MelderInfo_writeLine (U"Standard deviation in channel ", channel, U": ", Melder_single (stdev), U" Pascal");
		}
	}
}

double structSound :: v_getMatrix (const integer irow, const integer icol) const {
	if (irow < 1 || irow > our ny) {
		if (irow == 0) {
			if (icol < 1 || icol > nx)
				return 0.0;
			if (our ny == 1)
				return our z [1] [icol];   // optimization
			if (our ny == 2)
				return 0.5 * (our z [1] [icol] + our z [2] [icol]);   // optimization
			longdouble sum = 0.0;
			for (integer channel = 1; channel <= ny; channel ++)
				sum += our z [channel] [icol];
			return (double) sum / our ny;
		}
		return 0.0;
	}
	if (icol < 1 || icol > nx)
		return 0.0;
	return our z [irow] [icol];
}

double structSound :: v_getFunction2 (const double x, const double y) const {
	const integer channel = Melder_ifloor (y);
	if (channel < 0 || channel > our ny || y != (double) channel)
		return 0.0;
	return v_getFunction1 (channel, x);
}

autoSound Sound_create (integer numberOfChannels, double xmin, double xmax, integer nx, double dx, double x1) {
	try {
		autoSound me = Thing_new (Sound);
		Matrix_init (me.get(), xmin, xmax, nx, dx, x1, 1, numberOfChannels, numberOfChannels, 1, 1);
		return me;
	} catch (MelderError) {
		Melder_throw (U"Sound not created.");
	}
}

autoSound Sound_createSimple (integer numberOfChannels, double duration, double samplingFrequency) {
	Melder_assert (duration >= 0.0);
	Melder_assert (samplingFrequency > 0.0);
	const double numberOfSamples_real = round (duration * samplingFrequency);
	if (numberOfSamples_real > double (INT32_MAX))
		Melder_throw (U"Cannot create sounds with more than ", Melder_bigInteger (INT32_MAX), U" samples, because they cannot be saved to disk.");
	return Sound_create (numberOfChannels, 0.0, duration, (integer) (int32) numberOfSamples_real,
		1.0 / samplingFrequency, 0.5 / samplingFrequency);
}

autoSound Sound_convertToMono (constSound me) {
	if (my ny == 1)
		return Data_copy (me);   // optimization
	try {
		autoSound thee = Sound_create (1, my xmin, my xmax, my nx, my dx, my x1);
		if (my ny == 2) {   // optimization
			for (integer i = 1; i <= my nx; i ++)
				thy z [1] [i] = 0.5 * (my z [1] [i] + my z [2] [i]);
		} else {
			for (integer i = 1; i <= my nx; i ++) {
				longdouble sum = my z [1] [i] + my z [2] [i] + my z [3] [i];
				for (integer channel = 4; channel <= my ny; channel ++)
					sum += my z [channel] [i];
				thy z [1] [i] = double (sum / my ny);
			}
		}
		return thee;
	} catch (MelderError) {
		Melder_throw (me, U": not converted to mono.");
	}
}

autoSound Sound_convertToStereo (constSound me) {
	if (my ny == 2)
		return Data_copy (me);
	try {
		if (my ny > 2)
			Melder_throw (U"The Sound has ", my ny, U" channels; don't know which to choose.");
		Melder_assert (my ny == 1);
		autoSound thee = Sound_create (2, my xmin, my xmax, my nx, my dx, my x1);
		for (integer i = 1; i <= my nx; i ++)
			thy z [1] [i] = thy z [2] [i] = my z [1] [i];
		return thee;
	} catch (MelderError) {
		Melder_throw (me, U": not converted to stereo.");
	}
}

autoSound Sounds_combineToStereo (OrderedOf<structSound>* me) {
	try {
		integer totalNumberOfChannels = 0;
		double sharedSamplingPeriod = 0.0;
		for (integer isound = 1; isound <= my size; isound ++) {
			Sound sound = my at [isound];
			totalNumberOfChannels += sound -> ny;
			if (sharedSamplingPeriod == 0.0) {
				sharedSamplingPeriod = sound -> dx;
			} else if (sound -> dx != sharedSamplingPeriod) {
				Melder_throw (U"To combine sounds, their sampling frequencies should be equal.\n"
						U"You could resample one or more of the sounds before combining.");
			}
		}
		Melder_assert (my size > 0);
		double sharedMinimumTime = my at [1] -> xmin;
		double sharedMaximumTime = my at [1] -> xmax;
		for (integer isound = 2; isound <= my size; isound ++) {
			Sound sound = my at [isound];
			if (sound -> xmin < sharedMinimumTime)
				sharedMinimumTime = sound -> xmin;
			if (sound -> xmax > sharedMaximumTime)
				sharedMaximumTime = sound -> xmax;
		}
		autoINTVEC numberOfInitialZeroes = zero_INTVEC (my size);
		integer sharedNumberOfSamples = 0;
		double sumOfFirstTimes = 0.0;
		for (integer isound = 1; isound <= my size; isound ++) {
			Sound sound = my at [isound];
			numberOfInitialZeroes [isound] = Melder_ifloor ((sound -> xmin - sharedMinimumTime) / sharedSamplingPeriod);
			const double newFirstTime = sound -> x1 - sound -> dx * numberOfInitialZeroes [isound];
			sumOfFirstTimes += newFirstTime;
			const integer newNumberOfSamplesThroughLastNonzero = sound -> nx + numberOfInitialZeroes [isound];
			if (newNumberOfSamplesThroughLastNonzero > sharedNumberOfSamples)
				sharedNumberOfSamples = newNumberOfSamplesThroughLastNonzero;
		}
		const double sharedTimeOfFirstSample = sumOfFirstTimes / my size;   // this is an approximation
		autoSound thee = Sound_create (totalNumberOfChannels, sharedMinimumTime, sharedMaximumTime,
			sharedNumberOfSamples, sharedSamplingPeriod, sharedTimeOfFirstSample);
		integer channelNumber = 0;
		for (integer isound = 1; isound <= my size; isound ++) {
			Sound sound = my at [isound];
			const integer offset = numberOfInitialZeroes [isound];
			for (integer ichan = 1; ichan <= sound -> ny; ichan ++) {
				channelNumber ++;
				for (integer isamp = 1; isamp <= sound -> nx; isamp ++)
					thy z [channelNumber] [isamp + offset] = sound -> z [ichan] [isamp];
			}
		}
		return thee;
	} catch (MelderError) {
		Melder_throw (U"Sounds not combined to stereo.");
	}
}

static void Sound_checkChannelNumberWithinRange (constSound me, integer channel) {
	Melder_require (channel > 0,
		U"Your channel number is ", channel, U", but it should be positive.");
	Melder_require (channel <= my ny,
		U"Your channel number is ", channel,
		U", but it should not be greater than my number of channels, which is ",
		my ny, U"."
	);
}

autoSound Sound_extractChannel (constSound me, integer channel) {
	try {
		Sound_checkChannelNumberWithinRange (me, channel);
		autoSound you = Sound_create (1, my xmin, my xmax, my nx, my dx, my x1);
		your z.row (1)  <<=  my z.row (channel);
		return you;
	} catch (MelderError) {
		Melder_throw (me, U": channel ", channel, U" not extracted.");
	}
}

autoSound Sound_extractChannels (constSound me, constINTVECVU const& channels) {
	try {
		Melder_require (channels.size > 0,
			U"The number of channels should be greater than 0.");
		autoSound you = Sound_create (channels.size, my xmin, my xmax, my nx, my dx, my x1);
		for (integer ichan = 1; ichan <= channels.size; ichan ++) {
			Sound_checkChannelNumberWithinRange (me, channels [ichan]);
			your z.row (ichan)  <<=  my z.row (channels [ichan]);
		}
		return you;
	} catch (MelderError) {
		Melder_throw (me, U": channels not extracted.");
	}
}

static double getSumOfSquares (constSound me, double xmin, double xmax, integer *n) {
	Melder_assert (me);
	Function_unidirectionalAutowindow (me, & xmin, & xmax);
	integer imin, imax;
	*n = Sampled_getWindowSamples (me, xmin, xmax, & imin, & imax);
	if (*n <= 0)
		return undefined;
	longdouble sumOfSquares = 0.0;
	for (integer ichan = 1; ichan <= my ny; ichan ++) {
		constVECVU const& channel = my z.row (ichan);
		for (integer i = imin; i <= imax; i ++) {
			const longdouble value = channel [i];
			sumOfSquares += value * value;
		}
	}
	return double (sumOfSquares);
}

double Sound_getRootMeanSquare (constSound me, double xmin, double xmax) {
	integer n;
	const double sumOfSquares = getSumOfSquares (me, xmin, xmax, & n);
	return isdefined (sumOfSquares) ? sqrt (sumOfSquares / (n * my ny)) : undefined;
}

double Sound_getEnergy (constSound me, double xmin, double xmax) {
	integer n;
	const double sumOfSquares = getSumOfSquares (me, xmin, xmax, & n);
	return isdefined (sumOfSquares) ? sumOfSquares * my dx / my ny : undefined;
}

double Sound_getPower (constSound me, double xmin, double xmax) {
	integer n;
	const double sumOfSquares = getSumOfSquares (me, xmin, xmax, & n);
	return isdefined (sumOfSquares) ? sumOfSquares / (n * my ny) : undefined;
}

double Sound_getEnergyInAir (constSound me) {
	integer n;
	const double sumOfSquares = getSumOfSquares (me, 0.0, 0.0, & n);
	return isdefined (sumOfSquares) ? sumOfSquares * my dx / (400.0 * my ny) : undefined;
}

double Sound_getIntensity_dB (constSound me) {
	integer n;
	const double sumOfSquares = getSumOfSquares (me, 0.0, 0.0, & n);
	return isdefined (sumOfSquares) && sumOfSquares != 0.0 ? 10.0 * log10 (sumOfSquares / (n * my ny) / 4.0e-10) : undefined;
}

double Sound_getPowerInAir (constSound me) {
	integer n;
	const double sumOfSquares = getSumOfSquares (me, 0, 0, & n);
	return ( isdefined (sumOfSquares) ? sumOfSquares / (n * my ny) / 400.0 : undefined );
}

autoSound Matrix_to_Sound_mono (constMatrix me, integer rowNumber) {
	try {
		autoSound thee = Sound_create (1, my xmin, my xmax, my nx, my dx, my x1);
		if (rowNumber < 0)
			rowNumber = my ny + 1 + rowNumber;
		Melder_clip (1_integer, & rowNumber, my ny);
		thy z.row (1)  <<=  my z.row (rowNumber);
		return thee;
	} catch (MelderError) {
		Melder_throw (me, U": not converted to Sound.");
	}
}

autoSound Matrix_to_Sound (constMatrix me) {
	try {
		autoSound thee = Thing_new (Sound);
		my structMatrix :: v1_copy (thee.get());
		return thee;
	} catch (MelderError) {
		Melder_throw (me, U": not converted to Sound.");
	}
}

autoMatrix Sound_to_Matrix (constSound me) {
	try {
		autoMatrix thee = Thing_new (Matrix);
		my structMatrix :: v1_copy (thee.get());
		return thee;
	} catch (MelderError) {
		Melder_throw (me, U": not converted to Matrix.");
	}
}

autoSound Sound_upsample (constSound me) {
	try {
		constexpr integer antiTurnAround = 1000;
		constexpr integer sampleRateFactor = 2;
		constexpr integer numberOfPaddingSides = 2;   // namely beginning and end
		const integer nfft = Melder_iroundUpToPowerOfTwo (my nx + antiTurnAround * numberOfPaddingSides);
		const double newDx = my dx / sampleRateFactor;
		/*
			The computation of the new x1 relies on the idea that the left edge
			of the old first sample should coincide with the left edge of the new first sample
			(typically, e.g. if the old first sample starts at zero, which is usual,
			then the new first sample should also start at zero):
			old x1 - 0.5 * old dx == new x1 - 0.5 * new dx
			==>
			new x1 == old x1 - 0.5 * (old dx - new dx)
		*/
		autoSound thee = Sound_create (my ny, my xmin, my xmax, my nx * sampleRateFactor,
				newDx, my x1 - 0.5 * (my dx - newDx));
		for (integer ichan = 1; ichan <= my ny; ichan ++) {
			autoVEC data = zero_VEC (sampleRateFactor * nfft);   // zeroing is important...
			data.part (antiTurnAround + 1, antiTurnAround + my nx)  <<=  my z.row (ichan);   // ...because this fills only part of the sound
			NUMrealft (data.part (1, nfft), 1);
			const integer imin = (integer) (nfft * 0.95);
			for (integer i = imin + 1; i <= nfft; i ++)
				data [i] *= ((double) (nfft - i)) / (nfft - imin);
			data [2] = 0.0;
			NUMrealft (data.get(), -1);
			const double factor = 1.0 / nfft;
			for (integer i = 1; i <= thy nx; i ++)
				thy z [ichan] [i] = data [i + sampleRateFactor * antiTurnAround] * factor;
		}
		return thee;
	} catch (MelderError) {
		Melder_throw (me, U": not upsampled.");
	}
}

autoSound Sound_resample (constSound me, double samplingFrequency, integer precision) {
	const double upfactor = samplingFrequency * my dx;
	if (fabs (upfactor - 2.0) < 1e-6)
		return Sound_upsample (me);
	if (fabs (upfactor - 1.0) < 1e-6)
		return Data_copy (me);
	try {
		const integer numberOfSamples = Melder_iround ((my xmax - my xmin) * samplingFrequency);
		if (numberOfSamples < 1)
			Melder_throw (U"The resampled Sound would have no samples.");
		autoSound filtered;
		const bool weNeedAnAntiAliasingFilter = ( upfactor < 1.0 );
		if (weNeedAnAntiAliasingFilter) {
			constexpr integer antiTurnAround = 1000;
			constexpr integer numberOfPaddingSides = 2;   // namely beginning and end
			const integer nfft = Melder_iroundUpToPowerOfTwo (my nx + antiTurnAround * numberOfPaddingSides);
			autoVEC data = raw_VEC (nfft);   // will be zeroed in every turn of the loop
			filtered = Sound_create (my ny, my xmin, my xmax, my nx, my dx, my x1);
			for (integer ichan = 1; ichan <= my ny; ichan ++) {
				for (integer i = 1; i <= nfft; i ++)
					data [i] = 0.0;
				data.part (antiTurnAround + 1, antiTurnAround + my nx)  <<=  my z.row (ichan);
				NUMrealft (data.get(), 1);   // go to the frequency domain
				for (integer i = Melder_ifloor (upfactor * nfft); i <= nfft; i ++)
					data [i] = 0.0;   // filter away high frequencies
				data [2] = 0.0;
				NUMrealft (data.get(), -1);   // return to the time domain
				const double factor = 1.0 / nfft;
				const VEC to = filtered -> z.row (ichan);
				for (integer i = 1; i <= my nx; i ++)
					to [i] = data [i + antiTurnAround] * factor;
			}
			me = filtered.get();   // reference copy; remove at end
		}
		autoSound thee = Sound_create (my ny, my xmin, my xmax, numberOfSamples, 1.0 / samplingFrequency,
				0.5 * (my xmin + my xmax - (numberOfSamples - 1) / samplingFrequency));
		for (integer ichan = 1; ichan <= my ny; ichan ++) {
			if (precision <= 1) {
				for (integer i = 1; i <= numberOfSamples; i ++) {
					const double x = Sampled_indexToX (thee.get(), i);
					const double index = Sampled_xToIndex (me, x);
					const integer leftSample = Melder_ifloor (index);
					const double fraction = index - leftSample;
					thy z [ichan] [i] = ( leftSample < 1 || leftSample >= my nx ? 0.0 :
							(1 - fraction) * my z [ichan] [leftSample] + fraction * my z [ichan] [leftSample + 1] );
				}
			} else {
				for (integer i = 1; i <= numberOfSamples; i ++) {
					const double x = Sampled_indexToX (thee.get(), i);
					const double index = Sampled_xToIndex (me, x);
					thy z [ichan] [i] = NUM_interpolate_sinc (my z.row (ichan), index, precision);
				}
			}
		}
		return thee;
	} catch (MelderError) {
		Melder_throw (me, U": not resampled.");
	}
}

autoSound Sounds_append (constSound me, double silenceDuration, constSound thee) {
	try {
		const integer nx_silence = Melder_iround (silenceDuration / my dx), nx = my nx + nx_silence + thy nx;
		if (my ny != thy ny)
			Melder_throw (U"The numbers of channels are not equal (e.g. one is mono, the other stereo).");
		if (my dx != thy dx)
			Melder_throw (U"The sampling frequencies are not equal.");
		autoSound him = Sound_create (my ny, 0.0, nx * my dx, nx, my dx, 0.5 * my dx);
		his z.verticalBand (1, my nx)  <<=  my z.all();
		his z.verticalBand (my nx + nx_silence + 1, nx)  <<=  thy z.all();
		return him;
	} catch (MelderError) {
		Melder_throw (me, U" & ", thee, U": not appended.");
	}
}

autoSound Sounds_concatenate (SoundList list, double overlapTime) {
	try {
		integer sharedNumberOfChannels = 0, totalNumberOfSamples = 0;
		double sharedTimeStep = 0.0;
		for (integer i = 1; i <= list->size; i ++) {
			const Sound sound = list->at [i];
			if (sharedNumberOfChannels == 0) {
				sharedNumberOfChannels = sound -> ny;
			} else if (sound -> ny != sharedNumberOfChannels) {
				Melder_throw (U"To concatenate sounds, their numbers of channels (mono, stereo) should be equal.");
			}
			if (sharedTimeStep == 0.0) {
				sharedTimeStep = sound -> dx;
			} else if (sound -> dx != sharedTimeStep) {
				Melder_throw (U"To concatenate sounds, their sampling frequencies should be equal.\n"
						U"You could resample one or more of the sounds before concatenating.");
			}
			totalNumberOfSamples += sound -> nx;
		}
		const integer numberOfSmoothingSamples = Melder_iround (overlapTime / sharedTimeStep);
		totalNumberOfSamples -= numberOfSmoothingSamples * (list->size - 1);
		Melder_require (totalNumberOfSamples > 0,
			U"The joint overlap should not exceed the total duration of the sounds.\nChoose a shorter overlap time.");
		autoSound thee = Sound_create (sharedNumberOfChannels, 0.0, totalNumberOfSamples * sharedTimeStep,
				totalNumberOfSamples, sharedTimeStep, 0.5 * sharedTimeStep);
		autoVEC smoother;
		if (numberOfSmoothingSamples > 0) {
			smoother = raw_VEC (numberOfSmoothingSamples);
			const double factor = NUMpi / numberOfSmoothingSamples;
			for (integer i = 1; i <= numberOfSmoothingSamples; i ++)
				smoother [i] = 0.5 - 0.5 * cos (factor * (i - 0.5));
		}
		integer sampleOffset = 0;
		for (integer i = 1; i <= list->size; i ++) {
			const Sound sound = list->at [i];
			if (numberOfSmoothingSamples > 2 * sound -> nx)
				Melder_throw (U"At least one of the sounds is shorter than twice the overlap time.\nChoose a shorter overlap time.");
			const bool thisIsTheFirstSound = ( i == 1 );
			const bool thisIsTheLastSound = ( i == list->size );
			const bool weNeedSmoothingAtTheStartOfThisSound = ! thisIsTheFirstSound;
			const bool weNeedSmoothingAtTheEndOfThisSound = ! thisIsTheLastSound;
			const integer numberOfSmoothingSamplesAtTheStartOfThisSound = ( weNeedSmoothingAtTheStartOfThisSound ? numberOfSmoothingSamples : 0 );
			const integer numberOfSmoothingSamplesAtTheEndOfThisSound = ( weNeedSmoothingAtTheEndOfThisSound ? numberOfSmoothingSamples : 0 );
			for (integer channel = 1; channel <= sharedNumberOfChannels; channel ++) {
				for (integer j = 1, mySample = 1, thySample = mySample + sampleOffset;
					 j <= numberOfSmoothingSamplesAtTheStartOfThisSound;
					 j ++, mySample ++, thySample ++)
				{
					thy z [channel] [thySample] += sound -> z [channel] [mySample] * smoother [j];   // add
				}
				thy z.row (channel).part (sampleOffset + 1 + numberOfSmoothingSamplesAtTheStartOfThisSound,
						sampleOffset + sound -> nx - numberOfSmoothingSamplesAtTheEndOfThisSound)  <<=
						sound -> z.row (channel).part (1 + numberOfSmoothingSamplesAtTheStartOfThisSound,
						sound -> nx - numberOfSmoothingSamplesAtTheEndOfThisSound);
				for (integer j = 1, mySample = sound -> nx - numberOfSmoothingSamplesAtTheEndOfThisSound + 1, thySample = mySample + sampleOffset;
					 j <= numberOfSmoothingSamplesAtTheEndOfThisSound;
					 j ++, mySample ++, thySample ++)
				{
					thy z [channel] [thySample] = sound -> z [channel] [mySample] * smoother [numberOfSmoothingSamplesAtTheEndOfThisSound + 1 - j];   // replace (or add, which is the same since it's all zeroes to start with)
				}
			}
			sampleOffset += sound -> nx - numberOfSmoothingSamplesAtTheEndOfThisSound;
		}
		Melder_assert (sampleOffset == totalNumberOfSamples);
		return thee;
	} catch (MelderError) {
		Melder_throw (U"Sounds not concatenated.");
	}
}

autoSound Sounds_convolve (constSound me, constSound thee, kSounds_convolve_scaling scaling, kSounds_convolve_signalOutsideTimeDomain signalOutsideTimeDomain) {
	try {
		if (my ny > 1 && thy ny > 1 && my ny != thy ny)
			Melder_throw (U"The numbers of channels of the two sounds have to be equal or 1.");
		if (my dx != thy dx)
			Melder_throw (U"The sampling frequencies of the two sounds have to be equal.");
		const integer n1 = my nx, n2 = thy nx;
		const integer n3 = n1 + n2 - 1;
		const integer nfft = Melder_iroundUpToPowerOfTwo (n3);
		autoVEC data1 = raw_VEC (nfft);
		autoVEC data2 = raw_VEC (nfft);
		integer numberOfChannels = std::max (my ny, thy ny);
		autoSound him = Sound_create (numberOfChannels, my xmin + thy xmin, my xmax + thy xmax, n3, my dx, my x1 + thy x1);
		for (integer channel = 1; channel <= numberOfChannels; channel ++) {
			VEC a = my z.row (my ny == 1 ? 1 : channel);
			for (integer i = n1; i > 0; i --)
				data1 [i] = a [i];
			for (integer i = n1 + 1; i <= nfft; i ++)
				data1 [i] = 0.0;
			a = thy z.row (thy ny == 1 ? 1 : channel);
			for (integer i = n2; i > 0; i --)
				data2 [i] = a [i];
			for (integer i = n2 + 1; i <= nfft; i ++)
				data2 [i] = 0.0;
			NUMrealft (data1.get(), 1);
			NUMrealft (data2.get(), 1);
			data2 [1] *= data1 [1];
			data2 [2] *= data1 [2];
			for (integer i = 3; i <= nfft; i += 2) {
				const double temp = data1 [i] * data2 [i] - data1 [i + 1] * data2 [i + 1];
				data2 [i + 1] = data1 [i] * data2 [i + 1] + data1 [i + 1] * data2 [i];
				data2 [i] = temp;
			}
			NUMrealft (data2.get(), -1);
			a = his z.row (channel);
			for (integer i = 1; i <= n3; i ++)
				a [i] = data2 [i];
		}
		switch (signalOutsideTimeDomain) {
			case kSounds_convolve_signalOutsideTimeDomain::ZERO: {
				// do nothing
			} break;
			case kSounds_convolve_signalOutsideTimeDomain::SIMILAR: {
				for (integer channel = 1; channel <= numberOfChannels; channel ++) {
					const VEC a = his z.row (channel);
					const double edge = std::min (n1, n2);
					for (integer i = 1; i < edge; i ++) {
						const double factor = edge / i;
						a [i] *= factor;
						a [n3 + 1 - i] *= factor;
					}
				}
			} break;
			//case kSounds_convolve_signalOutsideTimeDomain_PERIODIC: {
				// do nothing
			//} break;
			default: Melder_fatal (U"Sounds_convolve: unimplemented outside-time-domain strategy ", (int) signalOutsideTimeDomain);
		}
		switch (scaling) {
			case kSounds_convolve_scaling::INTEGRAL: {
				Vector_multiplyByScalar (him.get(), my dx / nfft);
			} break;
			case kSounds_convolve_scaling::SUM: {
				Vector_multiplyByScalar (him.get(), 1.0 / nfft);
			} break;
			case kSounds_convolve_scaling::NORMALIZE: {
				double normalizationFactor = Matrix_getNorm (me) * Matrix_getNorm (thee);
				if (normalizationFactor != 0.0)
					Vector_multiplyByScalar (him.get(), 1.0 / nfft / normalizationFactor);
			} break;
			case kSounds_convolve_scaling::PEAK_099: {
				Vector_scale (him.get(), 0.99);
			} break;
			default: Melder_fatal (U"Sounds_convolve: unimplemented scaling ", (int) scaling);
		}
		return him;
	} catch (MelderError) {
		Melder_throw (me, U" & ", thee, U": not convolved.");
	}
}

autoSound Sounds_crossCorrelate (constSound me, constSound thee, kSounds_convolve_scaling scaling, kSounds_convolve_signalOutsideTimeDomain signalOutsideTimeDomain) {
	try {
		Melder_require (my ny == 1 || thy ny == 1 || my ny == thy ny,
			U"The numbers of channels of the two sounds have to be equal or 1.");
		Melder_require (my dx == thy dx,
			U"The sampling frequencies of the two sounds have to be equal.");
		const integer numberOfChannels = std::max (my ny, thy ny);
		const integer n1 = my nx, n2 = thy nx, n3 = n1 + n2 - 1;
		const integer nfft = Melder_iroundUpToPowerOfTwo (n3);
		autoVEC data1 = raw_VEC (nfft);
		autoVEC data2 = raw_VEC (nfft);
		const double my_xlast = my x1 + (n1 - 1) * my dx;
		autoSound him = Sound_create (numberOfChannels, thy xmin - my xmax, thy xmax - my xmin, n3, my dx, thy x1 - my_xlast);
		for (integer channel = 1; channel <= numberOfChannels; channel ++) {
			VEC a = my z.row (my ny == 1 ? 1 : channel);
			for (integer i = n1; i > 0; i --)
				data1 [i] = a [i];
			for (integer i = n1 + 1; i <= nfft; i ++)
				data1 [i] = 0.0;
			a = thy z.row (thy ny == 1 ? 1 : channel);
			for (integer i = n2; i > 0; i --)
				data2 [i] = a [i];
			for (integer i = n2 + 1; i <= nfft; i ++)
				data2 [i] = 0.0;
			NUMrealft (data1.get(), 1);
			NUMrealft (data2.get(), 1);
			data2 [1] *= data1 [1];
			data2 [2] *= data1 [2];
			for (integer i = 3; i <= nfft; i += 2) {
				const double temp = data1 [i] * data2 [i] + data1 [i + 1] * data2 [i + 1];   // reverse me by taking the conjugate of data1
				data2 [i + 1] = data1 [i] * data2 [i + 1] - data1 [i + 1] * data2 [i];   // reverse me by taking the conjugate of data1
				data2 [i] = temp;
			}
			NUMrealft (data2.get(), -1);
			a = his z.row (channel);
			for (integer i = 1; i < n1; i ++)
				a [i] = data2 [i + (nfft - (n1 - 1))];   // data for the first part ("negative lags") is at the end of data2
			for (integer i = 1; i <= n2; i ++)
				a [i + (n1 - 1)] = data2 [i];   // data for the second part ("positive lags") is at the beginning of data2
		}
		switch (signalOutsideTimeDomain) {
			case kSounds_convolve_signalOutsideTimeDomain::ZERO: {
				// do nothing
			} break;
			case kSounds_convolve_signalOutsideTimeDomain::SIMILAR: {
				for (integer channel = 1; channel <= numberOfChannels; channel ++) {
					double * const a = & his z [channel] [0];
					double const edge = std::min (n1, n2);
					for (integer i = 1; i < edge; i ++) {
						const double factor = edge / i;
						a [i] *= factor;
						a [n3 + 1 - i] *= factor;
					}
				}
			} break;
			//case kSounds_convolve_signalOutsideTimeDomain_PERIODIC: {
				// do nothing
			//} break;
			default: Melder_fatal (U"Sounds_crossCorrelate: unimplemented outside-time-domain strategy ", (int) signalOutsideTimeDomain);
		}
		switch (scaling) {
			case kSounds_convolve_scaling::INTEGRAL: {
				Vector_multiplyByScalar (him.get(), my dx / nfft);
			} break;
			case kSounds_convolve_scaling::SUM: {
				Vector_multiplyByScalar (him.get(), 1.0 / nfft);
			} break;
			case kSounds_convolve_scaling::NORMALIZE: {
				const double normalizationFactor = Matrix_getNorm (me) * Matrix_getNorm (thee);
				if (normalizationFactor != 0.0)
					Vector_multiplyByScalar (him.get(), 1.0 / nfft / normalizationFactor);
			} break;
			case kSounds_convolve_scaling::PEAK_099: {
				Vector_scale (him.get(), 0.99);
			} break;
			default: Melder_fatal (U"Sounds_crossCorrelate: unimplemented scaling ", (int) scaling);
		}
		return him;
	} catch (MelderError) {
		Melder_throw (me, U" & ", thee, U": not cross-correlated.");
	}
}

autoSound Sound_autoCorrelate (constSound me, kSounds_convolve_scaling scaling, kSounds_convolve_signalOutsideTimeDomain signalOutsideTimeDomain) {
	try {
		const integer numberOfChannels = my ny, n1 = my nx, n2 = n1 + n1 - 1;
		const integer nfft = Melder_iroundUpToPowerOfTwo (n2);
		autoVEC data = raw_VEC (nfft);
		const double my_xlast = my x1 + (n1 - 1) * my dx;
		autoSound thee = Sound_create (numberOfChannels, my xmin - my xmax, my xmax - my xmin, n2, my dx, my x1 - my_xlast);
		for (integer channel = 1; channel <= numberOfChannels; channel ++) {
			VEC a = my z.row (channel);
			for (integer i = n1; i > 0; i --)
				data [i] = a [i];
			for (integer i = n1 + 1; i <= nfft; i ++)
				data [i] = 0.0;
			NUMrealft (data.get(), 1);
			data [1] *= data [1];
			data [2] *= data [2];
			for (integer i = 3; i <= nfft; i += 2) {
				data [i] = data [i] * data [i] + data [i + 1] * data [i + 1];
				data [i + 1] = 0.0;   // reverse me by taking the conjugate of data1
			}
			NUMrealft (data.get(), -1);
			a = thy z.row (channel);
			for (integer i = 1; i < n1; i ++)
				a [i] = data [i + (nfft - (n1 - 1))];   // data for the first part ("negative lags") is at the end of data
			for (integer i = 1; i <= n1; i ++)
				a [i + (n1 - 1)] = data [i];   // data for the second part ("positive lags") is at the beginning of data
		}
		switch (signalOutsideTimeDomain) {
			case kSounds_convolve_signalOutsideTimeDomain::ZERO: {
				// do nothing
			} break;
			case kSounds_convolve_signalOutsideTimeDomain::SIMILAR: {
				for (integer channel = 1; channel <= numberOfChannels; channel ++) {
					const VEC a = thy z.row (channel);
					const double edge = n1;
					for (integer i = 1; i < edge; i ++) {
						const double factor = edge / i;
						a [i] *= factor;
						a [n2 + 1 - i] *= factor;
					}
				}
			} break;
			//case kSounds_convolve_signalOutsideTimeDomain_PERIODIC: {
				// do nothing
			//} break;
			default: Melder_fatal (U"Sounds_autoCorrelate: unimplemented outside-time-domain strategy ", (int) signalOutsideTimeDomain);
		}
		switch (scaling) {
			case kSounds_convolve_scaling::INTEGRAL: {
				Vector_multiplyByScalar (thee.get(), my dx / nfft);
			} break;
			case kSounds_convolve_scaling::SUM: {
				Vector_multiplyByScalar (thee.get(), 1.0 / nfft);
			} break;
			case kSounds_convolve_scaling::NORMALIZE: {
				const double normalizationFactor = Matrix_getNorm (me) * Matrix_getNorm (me);
				if (normalizationFactor != 0.0)
					Vector_multiplyByScalar (thee.get(), 1.0 / nfft / normalizationFactor);
			} break;
			case kSounds_convolve_scaling::PEAK_099: {
				Vector_scale (thee.get(), 0.99);
			} break;
			default: Melder_fatal (U"Sounds_autoCorrelate: unimplemented scaling ", (int) scaling);
		}
		return thee;
	} catch (MelderError) {
		Melder_throw (me, U": autocorrelation not computed.");
	}
}

void Sound_draw (constSound me, Graphics g,
	double tmin, double tmax, double minimum, double maximum, bool garnish, conststring32 method)
{
	const bool timesAreReversed = ( tmin > tmax );
	if (timesAreReversed)
		std::swap (tmin, tmax);
	Function_bidirectionalAutowindow (me, & tmin, & tmax);
	/*
		Domain expressed in sample numbers.
	*/
	integer ixmin, ixmax;
	const integer n = Matrix_getWindowSamplesX (me, tmin, tmax, & ixmin, & ixmax);
	if (n < 1)
		return;
	/*
		Automatic vertical range.
	*/
	if (minimum == maximum) {
		Matrix_getWindowExtrema (me, ixmin, ixmax, 1, my ny, & minimum, & maximum);
		if (minimum == maximum) {
			minimum -= 1.0;
			maximum += 1.0;
		}
	}
	/*
		Set coordinates for drawing.
	*/
	Graphics_setInner (g);
	for (integer channel = 1; channel <= my ny; channel ++) {
		Graphics_setWindow (g, timesAreReversed ? tmax : tmin, timesAreReversed ? tmin : tmax,
			minimum - (my ny - channel) * (maximum - minimum),
			maximum + (channel - 1) * (maximum - minimum)
		);
		if (str32str (method, U"bars") || str32str (method, U"Bars")) {
			for (integer ix = ixmin; ix <= ixmax; ix ++) {
				const double x = Sampled_indexToX (me, ix);
				const double y = std::min (my z [channel] [ix], maximum);
				const double left = std::max (x - 0.5 * my dx, tmin);
				const double right = std::min (x + 0.5 * my dx, tmax);
				Graphics_line (g, left, y, right, y);
				Graphics_line (g, left, y, left, minimum);
				Graphics_line (g, right, y, right, minimum);
			}
		} else if (str32str (method, U"poles") || str32str (method, U"Poles")) {
			for (integer ix = ixmin; ix <= ixmax; ix ++) {
				const double x = Sampled_indexToX (me, ix);
				Graphics_line (g, x, 0, x, my z [channel] [ix]);
			}
		} else if (str32str (method, U"speckles") || str32str (method, U"Speckles")) {
			for (integer ix = ixmin; ix <= ixmax; ix ++) {
				const double x = Sampled_indexToX (me, ix);
				Graphics_speckle (g, x, my z [channel] [ix]);
			}
		} else {
			/*
				The default: draw as a curve.
			*/
			Graphics_function (g, & my z [channel] [0], ixmin, ixmax,
					Matrix_columnToX (me, ixmin), Matrix_columnToX (me, ixmax));
		}
	}
	Graphics_setWindow (g, timesAreReversed ? tmax : tmin, timesAreReversed ? tmin : tmax, minimum, maximum);
	if (garnish && my ny == 2)
		Graphics_line (g, tmin, 0.5 * (minimum + maximum), tmax, 0.5 * (minimum + maximum));
	Graphics_unsetInner (g);
	if (garnish) {
		Graphics_drawInnerBox (g);
		Graphics_textBottom (g, true, U"Time (s)");
		Graphics_marksBottom (g, 2, true, true, false);
		Graphics_setWindow (g, tmin, tmax, minimum - (my ny - 1) * (maximum - minimum), maximum);
		Graphics_markLeft (g, minimum, true, true, false, nullptr);
		Graphics_markLeft (g, maximum, true, true, false, nullptr);
		if (minimum != 0.0 && maximum != 0.0 && (minimum > 0.0) != (maximum > 0.0))
			Graphics_markLeft (g, 0.0, true, true, true, nullptr);
		if (my ny == 2) {
			Graphics_setWindow (g, timesAreReversed ? tmax : tmin, timesAreReversed ? tmin : tmax,
					minimum, maximum + (my ny - 1) * (maximum - minimum));
			Graphics_markRight (g, minimum, true, true, false, nullptr);
			Graphics_markRight (g, maximum, true, true, false, nullptr);
			if (minimum != 0.0 && maximum != 0.0 && (minimum > 0.0) != (maximum > 0.0))
				Graphics_markRight (g, 0.0, true, true, true, nullptr);
		}
	}
}

static double interpolate (constSound me, integer i1, integer channel)
/* Precondition: my z [1] [i1] != my z [1] [i1 + 1]; */
{
	const integer i2 = i1 + 1;
	const double x1 = Sampled_indexToX (me, i1), x2 = Sampled_indexToX (me, i2);
	const double y1 = my z [channel] [i1], y2 = my z [channel] [i2];
	return x1 + (x2 - x1) * y1 / (y1 - y2);   // linear
}
double Sound_getNearestZeroCrossing (constSound me, double position, integer channel) {
	VEC amplitude = my z.row (channel);
	const integer leftSample = Sampled_xToLowIndex (me, position);
	const integer rightSample = leftSample + 1;
	/* Are we already at a zero crossing? */
	if (leftSample >= 1 && rightSample <= my nx &&
		(amplitude [leftSample] >= 0.0) !=
		(amplitude [rightSample] >= 0.0))
	{
		return interpolate (me, leftSample, channel);
	}
	integer ileft, iright;
	double leftZero, rightZero;

	/*
		Search to the left.
	*/
	if (leftSample > my nx)
		return undefined;
	for (ileft = leftSample - 1; ileft >= 1; ileft --)
		if ((amplitude [ileft] >= 0.0) != (amplitude [ileft + 1] >= 0.0)) {
			leftZero = interpolate (me, ileft, channel);
			break;
		}

	/*
		Search to the right.
	*/
	if (rightSample < 1)
		return undefined;
	for (iright = rightSample + 1; iright <= my nx; iright ++)
		if ((amplitude [iright] >= 0.0) != (amplitude [iright - 1] >= 0.0)) {
			rightZero = interpolate (me, iright - 1, channel);
			break;
		}

	if (ileft < 1 && iright > my nx)
		return undefined;
	return ileft < 1 ? rightZero : iright > my nx ? leftZero :
		position - leftZero < rightZero - position ? leftZero : rightZero;
}

void Sound_setZero (mutableSound me, double tmin_in, double tmax_in, bool roundTimesToNearestZeroCrossing) {
	Function_unidirectionalAutowindow (me, & tmin_in, & tmax_in);
	Function_intersectRangeWithDomain (me, & tmin_in, & tmax_in);
	for (integer channel = 1; channel <= my ny; channel ++) {
		/* mutable */ double tmin = tmin_in, tmax = tmax_in;
		if (roundTimesToNearestZeroCrossing) {
			if (tmin > my xmin)
				tmin = Sound_getNearestZeroCrossing (me, tmin_in, channel);
			if (tmax < my xmax)
				tmax = Sound_getNearestZeroCrossing (me, tmax_in, channel);
		}
		if (isundef (tmin))
			tmin = my xmin;
		if (isundef (tmax))
			tmax = my xmax;
		/* mutable */ integer imin, imax;
		Sampled_getWindowSamples (me, tmin, tmax, & imin, & imax);
		for (integer i = imin; i <= imax; i ++)
			my z [channel] [i] = 0.0;
	}
}

autoSound Sound_createAsPureTone (integer numberOfChannels, double startingTime, double endTime,
	double sampleRate, double frequency, double amplitude, double fadeInDuration, double fadeOutDuration)
{
	try {
		Melder_require (numberOfChannels >= 1,
			U"The number of channels should be at least 1.");
		const double numberOfSamples_real = round ((endTime - startingTime) * sampleRate);
		if (numberOfSamples_real > double (INT32_MAX))
			Melder_throw (U"Cannot create sounds with more than ", Melder_bigInteger (INT32_MAX), U" samples, because they cannot be saved to disk.");
		autoSound me = Sound_create (numberOfChannels, startingTime, endTime, (integer) (int32) numberOfSamples_real,
			1.0 / sampleRate, startingTime + 0.5 / sampleRate);
		for (integer isamp = 1; isamp <= my nx; isamp ++) {
			const double time = my x1 + (isamp - 1) * my dx;
			/* mutable */ double value = amplitude * sin (2.0 * NUMpi * frequency * time);
			const double timeFromStart = time - startingTime;
			if (timeFromStart < fadeInDuration)
				value *= 0.5 - 0.5 * cos (NUMpi * timeFromStart / fadeInDuration);
			const double timeFromEnd = endTime - time;
			if (timeFromEnd < fadeOutDuration)
				value *= 0.5 - 0.5 * cos (NUMpi * timeFromEnd / fadeOutDuration);
			for (integer ichan = 1; ichan <= my ny; ichan ++)
				my z [ichan] [isamp] = value;
		}
		return me;
	} catch (MelderError) {
		Melder_throw (U"Sound not created from tone complex.");
	}
}

autoSound Sound_createAsToneComplex (double startTime, double endTime, double samplingFrequency,
	int phase, double frequencyStep, double firstFrequency, double ceiling, integer numberOfComponents)
{
	try {
		if (frequencyStep == 0.0)
			Melder_throw (U"The frequency step should not be zero.");
		/*
			Translate default firstFrequency.
		*/
		if (firstFrequency <= 0.0)
			firstFrequency = frequencyStep;
		const double firstOmega = 2.0 * NUMpi * firstFrequency;
		/*
			Translate default ceiling.
		*/
		const double omegaStep = 2.0 * NUMpi * frequencyStep;
		const double nyquistFrequency = 0.5 * samplingFrequency;
		if (ceiling <= 0.0 || ceiling > nyquistFrequency)
			ceiling = nyquistFrequency;
		/*
			Translate number of components.
		*/
		const integer maximumNumberOfComponents = Melder_ifloor ((ceiling - firstFrequency) / frequencyStep) + 1;
		if (numberOfComponents <= 0 || numberOfComponents > maximumNumberOfComponents)
			numberOfComponents = maximumNumberOfComponents;
		if (numberOfComponents < 1)
			Melder_throw (U"There would be zero sine waves.");
		/*
			Generate the Sound.
		*/
		const double factor = 0.99 / numberOfComponents;
		autoSound me = Sound_create (1, startTime, endTime, Melder_iround ((endTime - startTime) * samplingFrequency),
			1.0 / samplingFrequency, startTime + 0.5 / samplingFrequency);
		for (integer isamp = 1; isamp <= my nx; isamp ++) {
			const double t = Sampled_indexToX (me.get(), isamp);
			const double omegaStepT = omegaStep * t;
			const double firstOmegaT = firstOmega * t;
			longdouble value = 0.0;
			if (phase == Sound_TONE_COMPLEX_SINE) {
				for (integer icomp = 1; icomp <= numberOfComponents; icomp ++)
					value += sin (firstOmegaT + (icomp - 1) * omegaStepT);
			} else {
				for (integer icomp = 1; icomp <= numberOfComponents; icomp ++)
					value += cos (firstOmegaT + (icomp - 1) * omegaStepT);
			}
			my z [1] [isamp] = double (value) * factor;
		}
		return me;
	} catch (MelderError) {
		Melder_throw (U"Sound not created as tone complex.");
	}
}

void Sound_multiplyByWindow (mutableSound me, kSound_windowShape windowShape) {
	for (integer ichan = 1; ichan <= my ny; ichan ++) {
		const integer n = my nx;
		VEC channel = my z.row (ichan);
		switch (windowShape) {
			case kSound_windowShape::RECTANGULAR: {
				;
			} break; case kSound_windowShape::TRIANGULAR: {   // "Bartlett"
				for (integer i = 1; i <= n; i ++) {
					const double phase = double (i - 0.5) / n;   // 0..1
					channel [i] *= 1.0 - fabs ((2.0 * phase - 1.0));
				}
			} break; case kSound_windowShape::PARABOLIC: {   // "Welch"
				for (integer i = 1; i <= n; i ++) {
					const double phase = double (i - 0.5) / n;
					channel [i] *= 1.0 - (2.0 * phase - 1.0) * (2.0 * phase - 1.0);
				}
			} break; case kSound_windowShape::HANNING: {
				for (integer i = 1; i <= n; i ++) {
					const double phase = double (i - 0.5) / n;
					channel [i] *= 0.5 * (1.0 - cos (2.0 * NUMpi * phase));
				}
			} break; case kSound_windowShape::HAMMING: {
				for (integer i = 1; i <= n; i ++) {
					const double phase = double (i - 0.5) / n;
					channel [i] *= 0.54 - 0.46 * cos (2.0 * NUMpi * phase);
				}
			} break; case kSound_windowShape::GAUSSIAN_1: {
				const double imid = 0.5 * double (n + 1), edge = exp (-3.0), onebyedge1 = 1.0 / (1.0 - edge);
				for (integer i = 1; i <= n; i ++) {
					const double phase = ((double) i - imid) / n;   // -0.5..+0.5
					channel [i] *= (exp (-12.0 * phase * phase) - edge) * onebyedge1;
				}
			} break; case kSound_windowShape::GAUSSIAN_2: {
				const double imid = 0.5 * double (n + 1), edge = exp (-12.0), onebyedge1 = 1.0 / (1.0 - edge);
				for (integer i = 1; i <= n; i ++) {
					const double phase = ((double) i - imid) / n;
					channel [i] *= (exp (-48.0 * phase * phase) - edge) * onebyedge1;
				}
			} break; case kSound_windowShape::GAUSSIAN_3: {
				const double imid = 0.5 * double (n + 1), edge = exp (-27.0), onebyedge1 = 1.0 / (1.0 - edge);
				for (integer i = 1; i <= n; i ++) {
					const double phase = ((double) i - imid) / n;
					channel [i] *= (exp (-108.0 * phase * phase) - edge) * onebyedge1;
				}
			} break; case kSound_windowShape::GAUSSIAN_4: {
				const double imid = 0.5 * double (n + 1), edge = exp (-48.0), onebyedge1 = 1.0 / (1.0 - edge);
				for (integer i = 1; i <= n; i ++) {
					const double phase = ((double) i - imid) / n;
					channel [i] *= (exp (-192.0 * phase * phase) - edge) * onebyedge1;
				}
			} break; case kSound_windowShape::GAUSSIAN_5: {
				const double imid = 0.5 * double (n + 1), edge = exp (-75.0), onebyedge1 = 1.0 / (1.0 - edge);
				for (integer i = 1; i <= n; i ++) {
					const double phase = ((double) i - imid) / n;
					channel [i] *= (exp (-300.0 * phase * phase) - edge) * onebyedge1;
				}
			} break; case kSound_windowShape::KAISER_1: {
				const double imid = 0.5 * double (n + 1);
				const double factor = 1.0 / NUMbessel_i0_f (2 * NUMpi);
				for (integer i = 1; i <= n; i ++) {
					const double phase = 2.0 * ((double) i - imid) / n;   // -1..+1
					const double root = 1.0 - phase * phase;
					channel [i] *= root <= 0.0 ? 0.0 : factor * NUMbessel_i0_f (2.0 * NUMpi * sqrt (root));
				}
			} break; case kSound_windowShape::KAISER_2: {
				const double imid = 0.5 * double (n + 1);
				const double factor = 1.0 / NUMbessel_i0_f (2 * NUMpi * NUMpi + 0.5);
				for (integer i = 1; i <= n; i ++) {
					const double phase = 2.0 * ((double) i - imid) / n;   // -1..+1
					const double root = 1.0 - phase * phase;
					channel [i] *= root <= 0.0 ? 0.0 : factor * NUMbessel_i0_f ((2.0 * NUMpi * NUMpi + 0.5) * sqrt (root));
				}
			} break; default: {
			}
		}
	}
}

void Sound_scaleIntensity (mutableSound me, double newAverageIntensity) {
	const double currentIntensity = Sound_getIntensity_dB (me);
	if (isundef (currentIntensity))
		return;
	const double factor = pow (10.0, (newAverageIntensity - currentIntensity) / 20.0);
	my z.all()  *=  factor;
}

void Sound_overrideSamplingFrequency (mutableSound me, double rate) {
	my dx = 1.0 / rate;
	my x1 = my xmin + 0.5 * my dx;
	my xmax = my xmin + my nx * my dx;
}

autoSound Sound_extractPart (constSound me, double tmin, double tmax, kSound_windowShape windowShape, double relativeWidth, bool preserveTimes) {
	try {
		/*
			We do not clip to the Sound's time domain.
			Any samples outside it are taken to be zero.
		*/

		Function_unidirectionalAutowindow (me, & tmin, & tmax);
		/*
			Allow window tails outside specified domain.
		*/
		if (relativeWidth != 1.0) {
			const double margin = 0.5 * (relativeWidth - 1) * (tmax - tmin);
			tmin -= margin;
			tmax += margin;
		}
		/*
			Determine index range. We use all the real or virtual samples that fit within [t1..t2].
		*/
		const integer itmin = 1 + Melder_iceiling ((tmin - my x1) / my dx);
		const integer itmax = 1 + Melder_ifloor   ((tmax - my x1) / my dx);
		Melder_require (itmax >= itmin,
			U"Extracted Sound would contain no samples.");
		/*
			Create sound, optionally shifted to [0..t2-t1].
		*/
		autoSound thee = Sound_create (my ny, tmin, tmax, itmax - itmin + 1, my dx, my x1 + (itmin - 1) * my dx);
		if (! preserveTimes) {
			thy xmin = 0.0;
			thy xmax -= tmin;
			thy x1 -= tmin;
		}
		/*
			Copy only *real* samples into the new sound.
			The *virtual* samples will remain at zero.
		*/
		for (integer ichan = 1; ichan <= my ny; ichan ++) {
			const integer itmin_clipped = Melder_clippedLeft (1_integer, itmin);
			const integer itmax_clipped = Melder_clippedRight (itmax, my nx);
			thy z.row (ichan).part (1 - itmin + itmin_clipped, 1 - itmin + itmax_clipped)
					<<=  my z.row (ichan).part (itmin_clipped, itmax_clipped);
		}
		/*
			Multiply by a window that extends throughout the target domain.
		*/
		Sound_multiplyByWindow (thee.get(), windowShape);
		return thee;
	} catch (MelderError) {
		Melder_throw (me, U": part not extracted.");
	}
}

autoSound Sound_extractPartForOverlap (constSound me, double tmin, double tmax, double overlap) {
	try {
		Function_unidirectionalAutowindow (me, & tmin, & tmax);
		if (overlap > 0.0) {
			const double margin = 0.5 * overlap;
			tmin -= margin;
			tmax += margin;
		}
		Melder_clipLeft (my xmin, & tmin);
		Melder_clipRight (& tmax, my xmax);
		/*
			Determine index range. We use all the real or virtual samples that fit within [t1..t2].
		*/
		const integer itmin = 1 + Melder_iceiling ((tmin - my x1) / my dx);
		const integer itmax = 1 + Melder_ifloor   ((tmax - my x1) / my dx);
		Melder_require (itmax >= itmin,
			U"Extracted Sound would contain no samples.");
		/*
			Create sound.
		*/
		autoSound thee = Sound_create (my ny, tmin, tmax, itmax - itmin + 1, my dx, my x1 + (itmin - 1) * my dx);
		thy xmin = 0.0;
		thy xmax -= tmin;
		thy x1 -= tmin;
		/*
			Copy only *real* samples into the new sound.
			The *virtual* samples will remain at zero.
		*/
		for (integer ichan = 1; ichan <= my ny; ichan ++) {
			const integer itmin_clipped = Melder_clippedLeft (1_integer, itmin);
			const integer itmax_clipped = Melder_clippedRight (itmax, my nx);
			thy z.row (ichan).part (1 - itmin + itmin_clipped, 1 - itmin + itmax_clipped)
					<<=  my z.row (ichan).part (itmin_clipped, itmax_clipped);
		}
		return thee;
	} catch (MelderError) {
		Melder_throw (me, U": part not extracted.");
	}
}

void Sound_filterWithFormants (mutableSound me, double tmin, double tmax,
	int numberOfFormants, double formant [], double bandwidth [])
{
	try {
		Function_unidirectionalAutowindow (me, & tmin, & tmax);
		integer itmin, itmax;
		integer n = Sampled_getWindowSamples (me, tmin, tmax, & itmin, & itmax);
		Melder_require (n >= 3,
			U"Sound too short.");
		for (integer ichan = 1; ichan <= my ny; ichan ++) {
			VECVU const& channel = my z.row (ichan);
			VECVU const& window = channel.part (itmin, itmax);
			VECdeemphasize_f_inplace (window, my dx, 50.0);
			for (int iformant = 1; iformant <= numberOfFormants; iformant ++)
				VECfilterSecondOrderSection_fb_inplace (window, my dx, formant [iformant], bandwidth [iformant]);
		}
		Matrix_scaleAbsoluteExtremum (me, 0.99);
	} catch (MelderError) {
		Melder_throw (me, U": not filtered.");
	}
}

autoSound Sound_filter_oneFormant (constSound me, double frequency, double bandwidth) {
	try {
		autoSound thee = Data_copy (me);
		Sound_filterWithOneFormantInplace (thee.get(), frequency, bandwidth);
		return thee;
	} catch (MelderError) {
		Melder_throw (me, U": not filtered (one formant).");
	}
}

void Sound_filterWithOneFormantInplace (mutableSound me, double frequency, double bandwidth) {
	for (integer ichan = 1; ichan <= my ny; ichan ++) {
		VEC channel = my z.row (ichan);
		VECfilterSecondOrderSection_fb_inplace (channel, my dx, frequency, bandwidth);
	}
	Matrix_scaleAbsoluteExtremum (me, 0.99);
}

double Sound_computeEmphasisFactor (constSound me, double cutoffFrequency) {
	return exp (- NUM2pi * cutoffFrequency * my dx);
}

void Sound_preEmphasize_inplace (mutableSound me, double cutoffFrequency) {
	const double emphasisFactor = Sound_computeEmphasisFactor (me, cutoffFrequency);
	if (emphasisFactor != 0.0)   // OPTIMIZE; will happen for cut-off frequencies above 119 times the sampling frequency
		for (integer channel = 1; channel <= my ny; channel ++) {
			VEC s = my z.row (channel);
			for (integer i = my nx; i >= 2; i --)
				s [i] -= emphasisFactor * s [i - 1];
		}
}

void Sound_deEmphasize_inplace (Sound me, double cutoffFrequency) {
	const double emphasisFactor = Sound_computeEmphasisFactor (me, cutoffFrequency);
	if (emphasisFactor != 0.0)   // OPTIMIZE; will happen for cut-off frequencies above 119 times the sampling frequency
		for (integer channel = 1; channel <= my ny; channel ++) {
			VEC s = my z.row (channel);
			for (integer i = 2; i <= my nx; i ++)
				s [i] += emphasisFactor * s [i - 1];
		}
}

autoSound Sound_filter_preemphasis (constSound me, double frequency) {
	try {
		autoSound thee = Data_copy (me);
		Sound_preEmphasize_inplace (thee.get(), frequency);
		Matrix_scaleAbsoluteExtremum (thee.get(), 0.99);
		return thee;
	} catch (MelderError) {
		Melder_throw (me, U": not filtered (pre-emphasis).");
	}
}

autoSound Sound_filter_deemphasis (constSound me, double frequency) {
	try {
		autoSound thee = Data_copy (me);
		Sound_deEmphasize_inplace (thee.get(), frequency);
		Matrix_scaleAbsoluteExtremum (thee.get(), 0.99);
		return thee;
	} catch (MelderError) {
		Melder_throw (me, U": not filtered (de-emphasis).");
	}
}

void Sound_reverse (mutableSound me, double tmin, double tmax) {
	Function_unidirectionalAutowindow (me, & tmin, & tmax);
	integer itmin, itmax;
	const integer n = Sampled_getWindowSamples (me, tmin, tmax, & itmin, & itmax) / 2;
	for (integer channel = 1; channel <= my ny; channel ++) {
		double *amp = & my z [channel] [0];
		for (integer i = 0; i < n; i ++)
			std::swap (amp [itmin + i], amp [itmax - i]);
	}
}

autoSound Sounds_crossCorrelate_short (constSound me, constSound thee, double tmin, double tmax, bool normalize) {
	try {
		if (my dx != thy dx)
			Melder_throw (U"Sampling frequencies are not equal.");
		if (my ny != thy ny)
			Melder_throw (U"Numbers of channels are not equal.");
		const double dt = my dx;
		double dphase = (thy x1 - my x1) / dt;
		dphase -= Melder_roundDown (dphase);   // a number between 0 and 1
		const integer i1 = Melder_iceiling (tmin / dt - dphase);   // index of first sample if sample at dphase has index 0
		const integer i2 = Melder_ifloor   (tmax / dt - dphase);   // index of last sample if sample at dphase has index 0
		const integer nt = i2 - i1 + 1;
		if (nt < 1)
			Melder_throw (U"Window too small.");
		const double t1 = (dphase + i1) * dt;
		autoSound him = Sound_create (1, tmin, tmax, nt, dt, t1);
		for (integer i = 1; i <= nt; i ++) {
			const integer di = i - 1 + i1;
			for (integer ime = 1; ime <= my nx; ime ++) {
				if (ime + di < 1)
					continue;
				if (ime + di > thy nx)
					break;
				for (integer channel = 1; channel <= my ny; channel ++)
					his z [1] [i] += my z [channel] [ime] * thy z [channel] [ime + di];
			}
		}
		if (normalize) {
			longdouble mypower = 0.0, thypower = 0.0;
			for (integer channel = 1; channel <= my ny; channel ++) {
				for (integer i = 1; i <= my nx; i ++) {
					const double value = my z [channel] [i];
					mypower += value * value;
				}
				for (integer i = 1; i <= thy nx; i ++) {
					const double value = thy z [channel] [i];
					thypower += value * value;
				}
			}
			if (mypower != 0.0 && thypower != 0.0) {
				double factor = 1.0 / (sqrt (double (mypower)) * sqrt (double (thypower)));
				for (integer i = 1; i <= nt; i ++) {
					his z [1] [i] *= factor;
				}
			}
		} else {
			const double factor = dt / my ny;
			for (integer i = 1; i <= nt; i ++) {
				his z [1] [i] *= factor;
			}
		}
		return him;
	} catch (MelderError) {
		Melder_throw (me, U": not cross-correlated.");
	}
}

Thing_implement (SoundList, Ordered, 0);

/* End of file Sound.cpp */
